Utiliser OpenStreetMap avec R

Mise en pratique

Auteur·rice·s

Affiliation

Louis Laurian

Centre pour l’analyse spatiale et la géovisualisation (UAR RIATE)

Ronan Ysebaert

Timothée Giraud

Objectifs

Ce document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM): géocoder une adresse, import de tuiles et de couches géographiques, import de points d’intérêt, routage et calcul d’itinéraires.

Il propose aussi des méthodes spatiales (carroyages, KDE, segments de lignes) et des visualisations associées pour restituer graphiquement la structure spatiale de semis de points ou de lignes.

En utilisant exclusivement des informations issues d’OSM, l’exercice consiste à identifier les secteurs accessibles à pieds dans le voisinage géographique d’un lieu défini par une adresse.

Reproduire l’analyse

Pour reproduire l’analyse, télécharger l’archive ZIP associée au dépôt https://github.com/riatecom/osm-elementr-2025-application, ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation.

Cette analyse est reproductible en n’importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d’OSRM ont été paramétrés pour ne pas le surcharger. L’exécution de ces blocs de code prennent du temps.

En conséquence, 6 jeux de données ont été préparés en amont à partir de l’adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d’itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.

cs <- read.csv("data/case_studies.csv", encoding = "UTF-8", sep = ",")
cs

        lib                                                              adr
1  Grenoble                       1 Place de la Gare, 38000 Grenoble, France
2  Clermont  40 Avenue de l'Union Soviétique, 63000 Clermont-Ferrand, France
3  Toulouse                      Gare SNCF Matabiau, 31500 Toulouse, France 
4 Compiegne                       Gare de Compiègne, 60200 Compiègne, France
5 Condorcet                       Front Populaire, 93210 Saint-Denis, France
6    Nantes                27 Boulevard de Stalingrad, 44000 Nantes, France 
7     Lille                              Lille Flandres, 59000 Lille, France

Nous présentons ici les résultats depuis la gare SNCF de Compiègne.

cs <- cs[cs$lib == "Compiegne",]
adr <- cs$adr
lib <- cs$lib

0. Packages

# 1. Packages qui interfacent des données OSM
#install.packages('maposm', repos = 'https://riatelab.r-universe.dev')
library(tidygeocoder) # géocodage
library(maptiles) # import de tuiles OSM (raster)
library(osmdata) # import de données OSM (vecteur)
library(maposm) # import de données OSM (couches géo)
library(osrm) # calcul d'itinéraires

# 2. Autres packages utilitaires
library(sf) # manipulation de données vectorielles
library(terra) # manipulation de données raster
library(mapsf) # cartographie thématique
library(maplegend) # légendes
library(spatstat) # analyse de semis de points
library(mapview) # cartographie interactive
library(stplanr) # segmentiser des lignes
library(lwgeom)

1. Géocodage

Le package tidygeocoder (Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.

pt <- geo(address = adr) |>
  st_as_sf(coords = c("long", "lat"), crs = 4326)

Passing 1 address to the Nominatim single address geocoder

Query completed in: 1 seconds

pt

Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 2.82352 ymin: 49.42206 xmax: 2.82352 ymax: 49.42206
Geodetic CRS:  WGS 84
# A tibble: 1 × 2
  address                                              geometry
* <chr>                                             <POINT [°]>
1 Gare de Compiègne, 60200 Compiègne, France (2.82352 49.42206)

2. Import de tuiles OSM

Le package maptiles (Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster. 20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.

Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.

emprise <- st_buffer(pt, 3000) |>
  st_bbox()

tiles <- get_tiles(x = emprise,
                   project = FALSE,
                   crop = TRUE,
                   cachedir = "cache",
                   zoom = 14)

pt_reproj <- st_transform(pt, crs = "EPSG:3857")

mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = paste0("Gare de ", lib))
mf_title("Tuiles OSM - Défaut")

# Autre exemple de fournisseur
tiles <- get_tiles(x = emprise,
                    project = FALSE,
                    crop = TRUE,
                    provider = "OpenStreetMap.HOT",
                    zoom = 14)

mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = paste0("Gare de ", lib))
mf_title("Tuiles OSM - OpenStreetmap.HOT")

3. Import de couches géographiques OSM

Pour la cartographie, le package maposm (pas sur le CRAN) se base sur les fonctions du package osmdata (Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.

La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d’un point, dans un rayon donné. Le résultat prend la forme d’une liste d’objets sf.

res <- om_get(x = c(st_coordinates(pt)[1],
                    st_coordinates(pt)[2]),
              r = 2000)

Getting urban areas: 0.898 sec elapsed
Getting buildings: 6.727 sec elapsed
Getting green areas: 1.053 sec elapsed
Getting roads: 0.411 sec elapsed
Getting streets: 1.259 sec elapsed
Getting railways: 0.355 sec elapsed
Getting water bodies: 1.8 sec elapsed

names(res)

[1] "zone"     "urban"    "building" "green"    "road"     "street"   "railway" 
[8] "water"   

mf_map(res$zone, col = "#f2efe9", border = NA, add = FALSE)
mf_map(res$urban, col = "#e0dfdf", border = "#e0dfdf", lwd = .5, add = TRUE)
mf_map(res$green, col = "#c8facc", border = "#c8facc", lwd = .5, add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$building, col = "#d9d0c9", border = "#c6bab1", lwd = .5, add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Autour de la gare")
mf_arrow()
mf_scale(size = 500, scale_units = "m")
mf_credits("OpenStreetMap contributors, 2025", pos = "bottomleft")

4. Import de tags avec osmdata

Le package osmdata permet d’extraire des données OSM à partir de l’API Overpass. Il nécessite dans un premier temps d’établir une bounding box au sein de laquelle seront extraites les données, puis de l’exécuter (fonction add_osm_feature()). Le résultat est une liste, qu’il faut ensuite convertir en objet sf.

Les objets OSM prennent la forme de nodes (points), ways (lignes ou polygones), ou relations (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de tags, recensées sous forme de clé-valeurs (key-values). Ici, nous recherchons dans la clé amenity les points (nodes) correspondant aux restaurants, bars et cafés (value = restaurant, bar, cafe). En dehors des clés-valeurs principales, d’autres attributs peuvent être présents, comme en l’occurence le type de cuisine, ou bien les horaires.

# Définition d'une bounding box (emprise Compiègne)
q <- opq(bbox = emprise, osm_types = "node")

# Extraction des restaurants, bars et cafés
req <- add_osm_feature(opq = q,
                       key = 'amenity',
                       value = c("bar", "cafe", "restaurant"))

poi <- osmdata_sf(req)
poi <- poi$osm_points

# Sélectionner les variables
poi <- poi[, c("osm_id", "name", "amenity", "cuisine")]

# Intersection
poi <- st_intersection(poi, st_transform(res$zone, crs = 4326))

# Reprojection
poi_reproj  <- st_transform(poi , crs = "EPSG:3857")

La fonction om_map() affiche les couches géographiques récupérées.

# Affichage
om_map(res, title = "Autour de la gare")
mf_map(poi_reproj, type = "typo", var = "amenity", cex = .7, pch = 15, 
       leg_pos = "topright", leg_val_cex = .9, leg_title_cex = .9, add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, lwd = 2, col = "darkblue", add = TRUE)

Il est possible d’utiliser le package mapview (Appelhans et al. 2023) pour afficher les points ainsi que leurs attributs, de manière interactive.

mapview(res$zone, alpha.regions = 0, layer.name = "Zone") + 
  mapview(poi, layer.name = "Manger et boire")

osmextract

Le package osmextract (Gilardi et Lovelace 2024) permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d’abord télécharger l’ensemble des données OSM de la Picardie, disponibles ici, et pesant plus de 240 Mb.

Le code ci-dessous effectue les mêmes manipulations que précédemment.

library(osmextract)
osm_pt <- oe_get(place = "Picardie",      
                 layer = "points", 
                 extra_tags = "amenity", 
                 quiet = TRUE)

osm_pt <- osm_pt[osm_pt$amenity %in% c("bar", "cafe", "restaurant"), ]
osm_pt <- st_transform(osm_pt, crs = "EPSG:3857")
osm_pt <- st_filter(osm_pt, res$zone, .predicates = st_intersects)

5. Carroyage

Une fois les points d’intérêt récupérés, il est possible de les agréger dans une maille.

# Création d'une grille régulière hexagonale de 300m de résolution
grid <- st_make_grid(res$zone, cellsize = 300, square = FALSE)
grid <- st_intersection(grid, res$zone)
grid <- st_sf(ID = 1:length(grid), geom = grid)

# Compter les POIs dans les carreaux
inter <- st_intersects(grid, poi_reproj, sparse = TRUE)
grid$n_poi <- lengths(inter)

# Cartographie
om_map(res)
mf_map(grid, var = "n_poi", type = "choro", breaks = c(0,1,3,5,10,15),
       alpha = .6, border = NA, leg_title = "Nombre \nd'aménités",
       leg_val_rnd = 0, leg_pos = "topright", add = TRUE)
mf_title("Aménités gustatives - données carroyées")

6. Lissage (KDE)

Le package spatstat (Baddeley, Turner, et Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe ppp (planar point pattern), dans des fenêtres d’observation, de classe owin (observation window).

La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).

Moraga (2023)

p <- as.ppp(X = st_coordinates(poi_reproj), W = as.owin(res$zone))
ds <- density.ppp(x = p, sigma = 100, eps = 10, positive = TRUE)

Une fois les densités calculées, nous créons un objet de type SpatRaster. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.

# Calcul densité de restaurants par hectare
r <- rast(ds) * 100 * 100
crs(r) <- st_crs(poi_reproj)$wkt
r

class       : SpatRaster 
dimensions  : 400, 400, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 312312.9, 316312.9, 6344781, 6348781  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
name        :         lyr.1 
min value   : 2.225074e-304 
max value   :  1.968529e+00 

pal <- c( "#418e40", "#a5be6a", "#e2cf6a", "#c19849", "#ab6c32", "#842f1a")
mf_map(res$zone, col = "#f2efe9", border = NA, expandBB = c(0,.3,0,0))
mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left", 
          pal = pal, add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("La montagne culinaire - Compiègne")

Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction as.polygons() arrondit les valeurs des pixels à l’entier.

# par défaut
r_poly <- as.polygons(r) |>
  st_as_sf()

# tous les dixièmes
r_poly <- as.polygons(r, round = TRUE, digits = 1) |>
  st_as_sf()

mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5,
       add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r_poly, var = "lyr.1", breaks = quantile(r_poly$lyr.1),
       border = NA, lwd = .2, type = "choro", alpha = 0.6,
       leg_title = "Densité\n (n/ha)", add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)

Par aménité, discrétisation commune

7. Temps d’accès

Le package osrm (Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. Ces calculs s’effectuent sur la base de profils .lua qui pour chaque mode de transport définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil voiture est disponible ici.

D’après la documentation de l’API OSRM, le nombre de requêtes sur le serveur de démo est limité à 512 par connection à l’API, espacées d’au moins 5 secondes. Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).

7.1. Depuis la gare

Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la gare de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.

Nous construisons d’abord une matrice O/D composées de coordonnées X/Y des points d’origine et de destination.

Utilisation de l’API OSRM

Par souci de respect d’utilisation de l’API, une partie du code ci-dessous est passé en commentaire et non joué. Le résultat est enregistré sous format RDS pour être utilisé par la suite. Dans les morceaux de code suivants, le code apparaît dans des chunks sous option eval = FALSE.

# Création de la matrice O/D
o <- data.frame(X = st_coordinates(pt)[, 1],
                Y = st_coordinates(pt)[, 2])
row.names(o) <- 1:nrow(o)

d <- data.frame(X = st_coordinates(poi)[, 1],
                Y = st_coordinates(poi)[, 2])

# 1 requête : 1 point d'origine, 92 points de destination
#mat_dist <- osrmTable(src = o,
#                      dst = d,
#                      measure = "distance",
#                      osrm.profile = 'foot')
#saveRDS(mat_dist, "data/mat_dist.rds")

# Intervalle de 5 secondes entre les requêtes
#Sys.sleep(5)
#mat_dur <- osrmTable(src = o,
#                     dst = d,
#                     measure = "duration",
#                     osrm.profile = 'foot')
#saveRDS(mat_dur, "data/mat_dur.rds")

mat_dist <- readRDS("data/mat_dist.rds")
mat_dur <- readRDS("data/mat_dur.rds")

dist <- t(mat_dist$distances)
dur <- t(mat_dur$durations)
poi$dist <- dist
poi$dur <- dur

apply(poi$dur, 2, summary)

                1
Min.     1.000000
1st Qu.  7.600000
Median  10.500000
Mean     9.679747
3rd Qu. 11.750000
Max.    17.600000

apply(poi$dist, 2, summary)

                1
Min.      73.0000
1st Qu.  570.5000
Median   789.0000
Mean     729.0759
3rd Qu.  881.5000
Max.    1501.0000

Il est possible d’atteindre 40 restaurants, cafés ou bars en moins de 11 min à pieds depuis la gare de Compiègne, et les trois quarts se situent à moins de 880m à pieds.

7.2. Depuis la grille

Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille.

nrow(grid) * nrow(d)

[1] 14931

Nous avons en tout 14931 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination.

Le code pour effectuer les calculs est présenté ci-dessous, mais a été joués en local. Les résultats sont enregistrés dans le dossier data.

osrmNearest()

Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).

# Coordonnées des points d'origine : centroïdes des cellules de la grille
o <- st_centroid(grid) |>
  st_transform(crs = 4326)

o$X <- st_coordinates(o)[, 1]
o$Y <- st_coordinates(o)[, 2]
o <- st_set_geometry(o[, c("X", "Y")], NULL)

# Les points de destination restent les mêmes
matfoot <- list()

for(i in 1:nrow(d)){
  matfoot[[i]] <- osrmTable(src = o[, c("X", "Y")],
                        dst = d[i , c("X", "Y")],
                        measure = "duration",
                        osrm.profile = 'foot')
  Sys.sleep(5)
  matfoot[[i]]$ID_ORIG <- c(1:nrow(o))
  matfoot[[i]]$ID_DEST <- rep(i, nrow(o))
  matfoot[[i]] <- data.frame("ID_ORIG" = matfoot[[i]]$ID_ORIG,
                         "ID_DEST" = matfoot[[i]]$ID_DEST,
                         "duration" = matfoot[[i]]$durations)
  colnames(matfoot[[i]]) <- c("ID_ORIG", "ID_DEST", "dur_foot")
}

matfoot <- do.call(rbind, matfoot)
write.csv(matfoot, "data/matfoot.csv", row.names = FALSE)

Nous pouvons maintenant calculer plusieurs indicateurs d’accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.

mat <- read.csv("data/matfoot.csv")

# Restaurant le plus proche pour chaque point de grille
poi_min <- aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm = TRUE)

# Récupérer l'identifiant de destination associé
poi_min <- do.call(rbind,
                   lapply(split(mat, mat$ID_ORIG),
                          function(x) x[which.min(x$dur_foot), ]))

# Cartographier
poi_min_sf <- merge(grid, poi_min,
                    by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)
par(mfrow = c(1,1))
om_map(res, title = "Temps d'accès au restaurant le plus proche")
mf_map(poi_min_sf, var = "dur_foot", type = "choro", breaks = c(0,5,10,15,20,40),
       alpha = .6, border = NA, leg_title = "Temps\n(min.)",
       leg_val_rnd = 0, pal = "Greens", leg_pos = "topright", add = TRUE)

Valhalla

Valhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de tuiles contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu’elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, Contraction Hierarchies pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.

# moins de 5 minutes
mat5 <- mat[mat$dur_foot < 5 ,]
mat5$n <- 1
mat5 <- aggregate(list(n_resto = mat5$n),
                  list(ID_ORIG = mat5$ID_ORIG),
                  FUN = sum)

grid <- merge(grid[, "ID"], mat5,
              by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)
grid$n_resto[is.na(grid$n_resto)] <- 0

om_map(res, title = "Restaurants à moins de 5 minutes")
mf_map(grid, var = "n_resto", type = "choro", breaks = c(0,1,5,10,30,max(grid$n_resto)),
       alpha = .6, border = NA, leg_title = "Nombre de\nrestaurants",
       leg_val_rnd = 0, pal = "Reds", leg_pos = "topright", add = TRUE)

Le package osrm propose aussi des isochrones, à l’aide de la fonction osrmIsochrone().

isos <- osrmIsochrone(loc = pt,
                      breaks = seq(0,45,5),
                      res = 60,
                      osrm.profile = "foot")

st_write(isos, "data/mat.gpkg", layer = "isos", delete_layer = TRUE)

isos <- st_read("data/mat.gpkg", layer = "isos", quiet = TRUE)
isos <- st_transform(isos, crs = "EPSG:3857")
isos <- st_intersection(isos, res$zone)

om_map(res, title = "Isochrones autour de la gare de Compiègne")
mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser",
      border = "white", leg_pos = NA, add = TRUE)

8. Itinéraires

Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.

routes <- list()
o <- data.frame(X = st_coordinates(pt)[, 1],
                Y = st_coordinates(pt)[, 2])
row.names(o) <- 1
for(i in 1:nrow(d)){
  routes[[i]] <- osrmRoute(src = o,
                           dst = d[i,],
                           overview = "full",
                           osrm.profile = 'foot')
   Sys.sleep(5)
}

routes <- do.call(rbind, routes)
st_write(routes, "data/mat.gpkg", layer = "routes", delete_layer = TRUE)

Le package stplanr (Lovelace, Ellison, et Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.

routes <- st_read("data/mat.gpkg", layer = "routes", quiet = TRUE)
routes <- st_transform(routes, crs = "EPSG:3857")
routes$n <- 1
routes$distance <- routes$distance * 1000

# La fonction overline permet de compter le nombre d'occurence des tronçons
seg <- overline(routes, attrib = "n", fun = sum)

routes <- routes[order(routes$distance, decreasing = TRUE), ]
l_seg <- list()

# La fonction line_segment découpe un objet sf LINESTRING en segments réguliers
for(i in 1:nrow(routes)){
  l_seg[[i]] <- line_segment(
    routes[i, ],
    segment_length = max(routes$distance) / 20)
}

pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")
pal2 <- paste0(pal, "70")

om_map(res, title = "Routes", theme = "dark")

for(i in nrow(routes):1){
  mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)
}

leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)),
    pal = pal, pos = "left", title = "Distance (m)")

om_map(res, title = "Routes", theme = "dark")
mf_map(seg, var = "n", type = "grad", col= "#00f5f5", leg_pos = NA, 
       breaks = "fisher", nbreaks = 5, lwd = c(1,3,6,10,15), add = TRUE)
mf_map(seg, var = "n", type = "grad", col= "#fff70c", breaks = "fisher", 
       nbreaks = 5, lwd = c(.5,1.5,2.5,6,10), leg_val_rnd = 0,
       leg_pos = "topright", leg_title = 'N. Routes', add = TRUE)

om_map(res, title = "Routes", theme = "dark")
mf_map(seg, var = "n", type = "grad", col= "#00f5f570", leg_pos = NA, 
       breaks = "fisher", nbreaks = 5, lwd = c(2,5,9,14,18), add = TRUE)
mf_map(seg, var = "n", type = "grad", col= "#00f5f5", breaks = "fisher", 
       nbreaks = 5, lwd = c(1,3,7,11,13), leg_val_rnd = 0,
       leg_pos = "topright", leg_title = 'N. Routes', add = TRUE)

Bonus : la tournée des pizzeria

Je fais une escale d’une demie-heure à Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps depuis la gare. Avec la fonction osrmTrip(), c’est possible !

Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.

# Coordonnées des pizzeria
pizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]

# Ajouter le point de départ : la gare
pizz <- rbind(pt[1, "geometry"], pizz)
pizz$id <- 1:nrow(pizz)

trip <- osrmTrip(loc = pizz,
                 overview = "full",
                 osrm.profile = "foot")

trip <- trip[[1]]$trip
st_write(trip, "data/mat.gpkg", layer = "trip", delete_layer = TRUE)

trip <- st_read("data/mat.gpkg", layer = "trip", quiet = TRUE)
trips <- st_transform(trip, crs = "EPSG:3857")

# Pour les labels : récupérer le début de la ligne
start <- lwgeom::st_startpoint(trips)
start <- do.call(rbind, start)
start <- data.frame(id = rownames(trips), X = start[,1], Y = start[,2])
start <- st_as_sf(start, coords = c("X", "Y"), crs = "EPSG:3857")

par(mfrow = c(1,1))
om_map(res, theme = "pizza", title = "La route de la pizza")
mf_map(trips, col = "black", lwd = 3, lty = 3, add = TRUE)
mf_map(start[start$id != 1, ], pch = 20, cex = 2, col = "black", add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_label(start, var = "id", cex = 1.2, overlap = TRUE, pos = 4, col = "black", halo = TRUE)

sum(trip$duration)

[1] 23.86833

D’après les données OSM, moins de 25 minutes suffisent pour faire la tournée des pizzas à Compiègne. On est large !!

Les références

Appelhans, Tim, Florian Detsch, Christoph Reudenbach, et Stefan Woellauer. 2023. Interactive Viewing of Spatial Data in R. https://CRAN.R-project.org/package=mapview.

Baddeley, Adrian, Rolf Turner, et Ege Rubak. 2025. spatstat: Spatial Point Pattern Analysis, Model-Fitting, Simulation, Tests. https://CRAN.R-project.org/package=spatstat.

Cambon, Jesse, Diego Hernangómez, Christopher Belanger, Daniel Possenriede, et Otto Hansen. 2025. tidygeocoder: An intuitive interface for getting data from geocoding services. https://CRAN.R-project.org/package=tidygeocoder.

Gilardi, Andrea, et Robin Lovelace. 2024. osmextract: Download and Import Open Street Map Data Extracts. https://CRAN.R-project.org/package=osmextract.

Giraud, Timothée. 2024. osrm: Interface Between R and the OpenStreetMap-Based Routing Service OSRM. https://CRAN.R-project.org/package=osrm.

———. 2025. maptiles: Download and Display Map Tiles. https://CRAN.R-project.org/package=maptiles.

Lovelace, Robin, Richard Ellison, et Malcolm Morgan. 2024. stplanr: Sustainable Transport Planning. https://CRAN.R-project.org/package=stplanr.

Moraga, Paula. 2023. « Spatial Statistics for Data Science: Theory and Practice with R. » https://www.paulamoraga.com/book-spatial/index.html.

Padgham, Mark, Bob Rudis, Robin Lovelace, Maëlle Salmon, et Joan Maspons. 2023. osmdata: Import ’OpenStreetMap’ Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.